Fruit Bromelain-Derived Peptide Potentially Restrains the Attachment of SARS-CoV-2 Variants to hACE2: A Pharmacoinformatics Approach

Before entering the cell, the SARS-CoV-2 spike glycoprotein receptor-binding domain (RBD) binds to the human angiotensin-converting enzyme 2 (hACE2) receptor. Hence, this RBD is a critical target for the development of antiviral agents. Recent studies have discovered that SARS-CoV-2 variants with mutations in the RBD have spread globally. The purpose of this in silico study was to determine the potential of a fruit bromelain-derived peptide. DYGAVNEVK. to inhibit the entry of various SARS-CoV-2 variants into human cells by targeting the hACE binding site within the RBD. Molecular docking analysis revealed that DYGAVNEVK interacts with several critical RBD binding residues responsible for the adhesion of the RBD to hACE2. Moreover, 100 ns MD simulations revealed stable interactions between DYGAVNEVK and RBD variants derived from the trajectory of root-mean-square deviation (RMSD), radius of gyration (Rg), and root-mean-square fluctuation (RMSF) analysis, as well as free binding energy calculations. Overall, our computational results indicate that DYGAVNEVK warrants further investigation as a candidate for preventing SARS-CoV-2 due to its interaction with the RBD of SARS-CoV-2 variants.


Introduction
As of 21 October 2021, the WHO (World Health Organization) had received reports of 241,886,635 confirmed coronavirus disease 2019 (COVID-19) cases worldwide, with 4,919,755 deaths. Although global vaccination campaigns are currently underway, it remains unclear how long the vaccine will provide immune defense against infection or if currently approved vaccinations will be sufficient to protect against emerging virus variants. Currently, numerous SARS-CoV-2 variants are emerging worldwide. The Centers for Disease Control and Prevention (CDC) have classified the SARS-CoV-2 variants as variants of interest (VOI), variants of concern (VOC), and variants of high consequence (VOHC). Additionally, the WHO has given Greek letters to SARS-CoV-2 variants.
Three variants that have risen to prominence in their respective countries and have sparked concern are the B.1.1.7, B.1.351, and P.1 lineages [1]. The B.1.1.7 lineage (alpha variant) was first described in the United Kingdom, while the B.1.351 lineage (beta variant) was initially reported in South Africa, and the P.1 lineage (gamma variant) was first reported in Brazil. Each of these three variants contains the N501Y mutation, which converts the amino acid asparagine (N) to tyrosine (Y) in the receptor-binding domain (RBD) subunit S1b of the glycoprotein spike. The B.1.1.7 lineage only has the N501Y mutation, while the B.1.351 and P.1 lineages have multiple mutations in the spike protein, including K417N, E484K, and N501Y [1][2][3][4][5]. Notably, SARS-CoV-2 mutations in the RBD were not confined to these three variants. To date, several variants with various mutations in different positions have circulated world-wide, including the California, New York, Scotland, Nigeria, and Indian variants [6][7][8][9][10].
To trigger cell entry and infection, the SARS-CoV-2 spike protein interacts directly with an enzyme known as the human angiotensin-converting enzyme 2 (hACE2) receptor. What is more concerning is that some of these variants have been shown to bind to the hACE2 receptor more effectively [2]. While hACE2 receptors are located on the surfaces of cells in a variety of tissues, they are particularly prevalent in the lungs [3,4]. Concerns have arisen since it is suspected that some of the COVID-19 vaccines currently in use are less effective against these variants [5][6][7]. As a result, researchers are also seeking COVID-19 antidotes.
Recently, several studies have been directed toward investigating a peptide (a small part of a protein) that can prevent the binding of the SARS-CoV-2 RBD to hACE2 [8][9][10]. For example, an antiviral peptide is a form of antiviral agent that is intended to be used as a therapeutic agent against a particular disease, for example, COVID-19. According to a report, 124 studies involving peptides were conducted in the search for an antidote for COVID-19. Of these, there were several clinical trials in the management of COVID-19, including immunomodulatory (5 trials), homeostasis regaining (8 trials), vaccination (9 trials), and antiviral activity trials (4 trials) [11].
Since peptides are more effective and precise than small-molecule drugs, they may be better tolerated [12]. Moreover, unlike other antiviral drugs, peptides have no toxicity to human cells [13,14]. However, despite their high therapeutic potential, some peptides have failed to reach clinical trials because of their toxicity (hemolytic activity) [15]. Peptide synthesis can also be easily implemented and tuned. On the other hand, small molecules often necessitate the creation of new time-consuming and expensive synthetic techniques [16]. Antiviral peptides outperform traditional antiviral drugs [17] because they are more effective, have a smaller molecular weight, and have fewer side effects [18]. To date, the FDA has approved more than 60 peptide-based drugs. However, more than 150 peptides are still undergoing advanced clinical trials [19]. Enfuvirtide is an HIV-1 fusion inhibitor linear synthetic peptide with a length of 36-amino acids and is an FDA-approved antiviral peptide [20]. Thus, peptides are molecules that can be tested against SARS-CoV-2 to potentially develop new drugs to treat COVID-19 [21].
Bioactive peptides are defined as peptide sequences contained within a protein that have a beneficial effect on bodily functions and/or have a beneficial effect on human health in addition to the protein's known nutritional value. They typically have a length of 3-20 amino acid residues [22]. Bioactive peptides can be obtained from a variety of sources, including animals, plants, and microorganisms. They can also be derived from a variety of proteins. When protein is consumed, it is digested in the digestive tract, resulting in the formation of peptides with numerous beneficial properties for the body. One of the proteins that has been widely researched and is known to have many health benefits is bromelain. It is a mixture of various cysteine proteinases with similar amino acid sequences and is found in pineapple fruits and stems [1,[23][24][25]. Bromelain has been shown to reduce the expression of ACE-2 and TMPRSS2 in VeroE6 cells, as well as to significantly reduce the expression of the S-Ectodomain of SARS-CoV-2 [26]. According to our previous in silico research, bromelain has a high binding affinity for various RBD variants and binds directly at the binding site between RBDs and hACE2. This suggests that bromelain has the potential to inhibit SARS-CoV-2 attachment to hACE2 [27]. In one study, the bromelain peptide biomarker DYGAVNEVK was found in the plasma of bromelain-treated mice. This peptide is one of bromelain's four proteolytically active proteins, which contribute to its therapeutic properties [28]. In the present study, we performed molecular docking analysis and an MD simulation study of the peptide DYGAVNEVK against several variants of the SARS-CoV-2 RBD.

Multiple Sequence Alignment of Wild-Type RBD and Its Variants
Multiple sequence alignment (MSA) techniques are a collection of algorithmic solutions for aligning evolutionarily related sequences. MSA of the amino acid sequences of the wild-type (WT) receptor-binding domain (RBD) and its variants was performed using the UCSF Chimera package (release 1.15) [29].

Three-Dimensional Structures of the Bromelain-Derived Peptide and the RBD Variants
The three-dimensional (3D) structures of the bromelain-derived peptide DYGAVNEVK and RBD variants were modeled and minimized using the SWISS-MODEL web server (https://swissmodel.expasy.org/interactive; accessed on 9 September 2021) [30,31]. Moreover, the two-dimensional (2D) structure was generated in PDBsum (http://www.ebi. ac.uk/thornton-srv/databases/pdbsum; accessed on 10 September 2021) [32]. The WT RBD was retrieved from a protein data bank with PDB ID 6M0J (https://www.rcsb.org/ structure/6M0J; accessed on 12 September 2021). The structure control of the variants' model structures was checked by MolProbity Structure Assessment and a Ramachandran plot.

Equilibrium Dissociation Constant Analysis
The equilibrium dissociation constant (K D ) is defined by the ratio of ligand receptor unbinding to binding rates. For each complex, the K D was predicted by PRODIGY (PROtein binDIng enerGY) (https://bianca.science.uu.nl/prodigy/; accessed on 16 September 2021) [38].

Analysis of MM-GBSA Free Energy
Using the molecular mechanics-generalized born surface area (MM-GBSA) approach, the binding-free energy (BFE) of the protein-peptide complex was calculated as the difference between the energy of the bound complex and that of the unbound protein and peptide. The calculations were performed using the HawkDock web server (http: //cadd.zju.edu.cn/hawkdock/; accessed on 17 September 2021) [39].

Analysis of the Complex Interface
The general content of the interface area resulting from molecular docking tests was analyzed using the PDBsum database for the structural analysis of 3D structures (EMBL-EBI; http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html; accessed on 18 September 2021). The 2D visualization of this interaction was generated using LigPlot+ [40].

Molecular Dynamics Simulations Study
MD simulations were performed using Gromacs version 2019.2 via WebGRO for macromolecular simulations (https://simlab.uams.edu/; accessed on 21 July 2021) [41]. The topology of the bromalin peptide and RBD protein complex was established by choosing the am-ber99sb-ildn [42] force field and the simple point-charge (SPC) water model. For the MD system, the triclinic water box was preferred and neutralized by adding the appropriate 0.15 M NaCl salt. The energy minimization of the created system was carried out in 5000 steps with the steepest descent integrator. The canonical ensemble NVT (moles (N), volume (V), and temperature (T)) equilibrium phase of the system was carried out at 300 K using the 0.3 ns duration V-rescale method [43], and the isothermal-isobaric ensemble NPT (moles (N), pressure (P), and temperature (T)) equilibrium phase was carried out using the Parrinello-Rahman method [44] at 0.3 ns under 1 atm of pressure. The MD simulation was performed using a leapfrog MD integrator to form 5000 frames with a duration of 100 ns. Root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), and radius of gyration (Rg) trajectory analyses were performed.

Binding-Free Energies Calculation
The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) BFE calculation is widely used to analyze the stability and bonding strength of protein-ligand, proteinpeptide, and protein-protein complexes [45]. The BFE consists of polar solvation energies, solvent accessible surface area (SASA) energy, and electrostatic and van der Waals interactions. In this study, binding free energy calculations for the RBD S1b units and bromelain peptide complexes of SARS-CoV-2 variants were performed using the MM-PBSA method with 60 frames between 40 and 70 ns from the MD trajectory. The 'MmPbSaStat python' script provided in g_mmpbsa was used for the average binding energy calculations [45,46].

Results and Discussion
The spread of SARS-CoV-2 variants across several continents is a significant source of concern for global human health. The variants are rapidly transmissible and quickly become prevalent in populations. Notably, the spike (S) protein has accumulated a large number of mutations, particularly within the amino-terminal domain (NTD) and the RBD. The emergence of mutations in this spike has direct implications for the high rate of viral infection since a conformational change in the RBD has resulted in stronger binding to the ACE2 receptor. Single amino acid substitutions should be monitored because they can have phenotypic consequences [47,48].
The MSA results for the RBD wild type and its variants are presented in Figure 1 and Table 1. These results show the positions of mutations in the RBD of SARS-CoV-2.
The following mutations were found in the RBD and are listed in Table 1: K417N (lysine, positive polar to asparagine, neutral polar); K417T (lysine, neutral polar to threonine, neutral polar); N439K (asparagine, neutral polar to lysine, positive polar); L452R (leucine, neutral nonpolar to arginine, positive polar); S477G (serine, neutral polar to glycine, neutral nonpolar); S477N (serine, neutral polar to asparagine, neutral polar), E484K (glutamate, negative polar to lysine, positive polar); E484Q (glutamate, negative polar to glutamine, neutral polar); and N501Y (asparagine, neutral polar to tyrosine, neutral polar). According to this explanation, since the K417T, S477N, and N501Y mutations do not change the charge or polarity of the amino acid residues, they hypothetically do not cause a conformational change in the RBD. However, since other mutations undergo changes in the charge and polarity of the amino acid residues, these are likely to change the conformation and have an impact on the attachment of the RBD to hACE2. If the mutation includes negatively charged, polar, and hydrophilic amino acids, there will be an increase in RBD stability [49].

The 3D and 2D Structures of Bromelain-Derived Peptide
The 3D and 2D structures of the bromelain-derived peptide are presented in Figure 2. The structure was generated from the sequence DYGAVNEVK (ASP-TYR-GLY-ALA-VAL-ASN-GLU-VAL-LYS). The 3D structure was used for peptide-protein docking analysis. The position of the sequence was demonstrated in the previous study [28].

Physicochemical Properties
The peptide has a molecular mass of 993.4751 Da. Ideally, small-molecule drugs with typical molecular weights of 500 Da or less are preferred for oral bioavailability [50]. However, if the peptide is between 500 and 1000 Da, there is still hope for further exploration, as it represents an enormous opportunity for those willing to explore new frontiers of drug discovery [51].
The peptide has an isoelectric point of 4. This means at a pH 4, the peptide carries no net electrical charge (that is, it is electrically neutral) based on the statistical mean. The overall net charge of the peptide at pH 7 is −1. The hydrophobicity scale for the peptide is 18.84 kcal/mol based on the free energy of transfer. The more positive the hydrophobicity scale, the more hydrophobic the molecule. Peptides containing a high proportion of hydrophobic amino acids will demonstrate a detrimental effect on their solubility in water. When designing soluble peptides, a good rule of thumb is to charge 1 out of every 5 amino acids. If this cannot be accomplished, non-essential amino acids in the peptide sequence can be replaced with charged residues [52]. This, of course, has the potential to alter the peptide's nature [53]. Consequently, substitutions should be carefully considered.
The amino acid composition and the number of peptide bonds present in a peptide are used to calculate the predictability of its molar extinction coefficient. The predicted molar extinction coefficient of the peptide is 1490 M −1 cm −1 . The N-terminal of the sequence considered is D (Asp). Based on this, the estimated half-life is 1.1 h in mammalian reticulocytes in an in vitro study. The instability index (II) is computed to be −13.18; therefore, this classifies the peptide as stable. Additionally, the peptide is water soluble. Good oral bioavailability will occur if a drug candidate exhibits optimal permeability and solubility at the absorption site [54].

Allergenicity and Toxicity Prediction
Allergic reactions are difficult to predict because they entail complicated interactions between a chemical (allergen) and the immune system [55]. The allergen initiates a Th2 response, which causes B cells to generate IgE and activates eosinophils [56]. Eosinophil accumulation in tissues can be extremely detrimental, as it results in inflammation and tissue damage [57]. According to the prediction, the studied peptide has the potential to be an allergen with a Tanimoto similarity index of 0.71. The higher the value, the closer the two structures are [34]. Additionally, the prediction indicates that the propensity for in vitro aggregation is 0.00.
The toxicity analysis considered several parameters, including the LD 50 , predicted toxicity class, hepatotoxicity, carcinogenicity, immunotoxicity, mutagenicity, and cytotoxicity. Toxicity levels are classified as follows: classes 1 and 2 (fatal if swallowed), class 3 (toxic if swallowed), class 4 (harmful if swallowed), class 5 (possibly harmful if swallowed), and class 6 (nontoxic) [58]. The predicated LD 50 of the peptide is 500 mg/kg, which classified the molecule in class 4, indicating it is harmful if swallowed. It is predicted to be non-hepatotoxic with a probability of 0.93. Additionally, it is not shown to be carcinogenic, immunotoxic, mutagenic, and cytotoxic, with probabilities of 0.64, 0.99, 0.79, and 0.69, respectively.

The Predicted IC 50 Value of Bromelain-Derived Peptide
The predicted IC50 value of the bromelain-derived peptide using the AVP-IC 50Pred server IC 50 was 40.67 µM (moderately effective, predicted by support vector machines (SVMs)) and 6.85 µM (effective, predicted by random forest (RF)). The IC 50 value is a quantitative measure of the amount of a molecule or drug required to inhibit up to half (50%) of a specific biological process. Since these values are only estimations, further in vitro research is required to validate them. In comparison, the AHB1 and AHB2 peptides have been shown to neutralize SARS-CoV-2 with IC 50 values of 35 and 16 nM, respectively [16]. The peptides P9R and PR had IC 50 values of 0.9 µg/mL and 2.4 µg/mL against SARS-CoV-2, respectively [59]. The amino acid composition of peptides may affect their inhibitory activity [60,61].

Analysis of the Interaction between Bromelain-Derived Peptide and RBD Variants
The HADDOCK scoring function combines various energies and buried surface area to arrive at a single numerical score. The HADDOCK score is defined as: 1.0 E vdw + 0.2 E elec + 1.0 E desol + 0.1 E AIR , where E vdw is an intermolecular van der Waals energy (kcal/mol), E elec is an intermolecular electrostatic energy (kcal/mol), E desol is a desolvation energy (kcal/mol), and E AIR is a restraint violation energy (kcal/mol). The scores of the complex predicted using the HADDOCK2.2 web server are listed in Table 2. The HADDOCK score determined for each interaction of the peptide bromelain and RBD after docking can be described as follows: and NG (−82.1 ± 6.0). Since lower HADDOCK scores show a higher affinity between the peptide and protein, the interaction formed will be stronger and more stable [62]. Hence, the interaction between bromelain-derived peptide and the NG variant shows the highest affinity, followed by the SA, SC, BR, and US variants.
Binding affinity is used to assess and rank the strength of the interactions formed, which is also calculated by the equilibrium dissociation constant. Thus, the lower the K D value, the higher the affinity [63]. The lowest K D value was obtained from the interaction of the bromelain-derived peptide with the NG variant (9. The 2D visualization of the interaction generated using LigPlot+ is presented in Figures 3 and 4, and Supplementary Materials (Figures S1-S8). Hydrogen bonds (H-bonds) are shown as green dotted lines. H-bonds play a critical role in drug-receptor interactions and in the structural integrity of a large number of biological molecules [64]. In addition to the van der Waals interaction in a complex, intermolecular hydrogen bonds contribute to the scoring function used to assess docking effectiveness [65]. Rathod et al. [8] discovered that the majority of studied peptides have a higher affinity for ACE2 than for the RBD residue binding motif. However, another study revealed that α-helical peptides extracted from the protease domain (PD) of ACE2 bind very specifically to SARS-CoV-2 and are stable, which implies that they can block the virus [9]. As a result, it is proposed that short peptides can be administered directly via inhalation to critical organs for SARS-CoV-2 infection, which offers an appealing alternative to traditional drug development [10]. The interaction studies were primarily concerned with the efficient binding of bromelainderived peptide with the receptor-binding motifs (RBMs) of RBDs. The results show that bromelain-derived peptide formed H-bonds with the active site residues of RBM, which suggests a good propensity for efficient binding to RBDs and the ability to inhibit virus attachment to the hACE2 receptor. The LigPlot+ analysis of bromelain-derived peptide and the RBD reveals the interaction of hydrogen bonds with residues on the active side of RBD at Gln493   According to the results of amino acid interactions, the bromelain-derived peptide was able to interact with the receptor-binding motif (RBM) of RBD by blocking unique residues designated as important in the binding of the human angiotensin-converting enzyme 2 (ACE2) cell receptor. The residues that are crucial for SARS-CoV-2 RBD binding to hACE2 are Gly446, Tyr449, Leu455, Phe486, Tyr491, Gln493, Gly496, Gln498, Thr500, Asn501, and Gly502 [66][67][68][69], as well as a salt bridge contributed by Lys417 [69]. The residues Phe486, Gln493, and Asn501 are the most important residues in the RBD identified by the hACE2 receptor on infected human cells because they facilitate RBD-hACE2 interaction [17,70]. Furthermore, it was reported that the amino acid substitutions S477G and S477N enhance the binding of the SARS-CoV-2 spike to the hACE2 receptor [71]. Table 3 summarizes the interaction between the amino acid residues of RBDs and bromelain-derived peptide, as determined by PDBsum analysis. With regards to the 3 most critical residues, i.e., Phe486, Gln493, and Asn501, the bromelain-derived peptide interacts with RBD WT via an H-bond on Gln493 and a hydrophobic interaction on Asn501. In particular, the BR and SA variants showed hydrophobic interactions with Gln493, while the SG and IN variants only showed hydrophobic interactions with Phe486. On the other hand, the US and NG variants contain H-bonds with Gln493 and hydrophobic interactions with Asn501. There is an H-bond on Gln493 and a hydrophobic interaction on Phe486 in the SN and SC variants. For Phe486 and Gln493, the UK variant demonstrates only hydrophobic interactions. Based on these findings, 8 of the 10 examined variants showed interactions with Gln493, including 5 variants with H-bonds and 3 variants with hydrophobic interactions. While 5 of the 10 variants interacted with Phe486 via hydrophobic interactions, only 2 variants interacted with Asn501 via hydrophobic interactions. Since hydrophobic interactions are entropy-driven at room temperature, they play a significant role in the docked complex binding affinity in a given solvent system [72]. The lower frequency of H-bonds is due to the binding pocket being more hydrophobic [73]. The peptide and RBD of the spike protein have a high number of H-bonds and hydrophobic interactions, thus indicating a strong interaction [21]. When compared to bromelain, the binding position of bromelain-derived peptide is generally similar, namely at the hACE2-RBD binding site. This is demonstrated by the similarity of the interacting amino acids between RBD WT and variants and these two molecules, particularly the critical residues Phe486, Gln493, and Asn501 [27]. Table 3. List of interacting amino acids between RBDs and bromelain-derived peptide. The position of the interacting residues in pocket and mutation sites are indicated in italics, while the key amino acid residues that play a role in binding RBD to hACE2 are marked in bold.
Proteins from mealworms, silkworm cocoons, and housefly larvae produce peptides when digested in the digestive tract. An in silico study reported that these peptides can bind to the RBD at key residues [75]. The peptide VEDKGMMHQQRMMEKAMNIPRM-CGTMQRKCRMS, derived from quinoa seed protein, was shown to form hydrophobic interactions with the key binding residues Leu455, Phe456, Phe486, and Gln493 on the RBD [76]. However, not all peptides can bind to key RBD residues. The chimeric peptides cnCoVP-2, cnCoVP-5, and cnCoVP-6 have been shown to interact with RBD, but not with the key residue (Phe486) [17]. Importantly, despite the substitution of neutral and polar asparagine for neutral and polar tyrosine in residue 501 in this study, the bromelain-derived peptide could still recognize and bind to it. As a result, the bromelain-derived peptide could potentially prevent the interaction between hACE2 and RBD that undergoes mutation at residue 501.

Prediction of the Position of Bromelain-Derived Peptide Inhibition between the RBD and hACE2
As illustrated in Figure 5, the type of bromelain-derived peptide inhibition at the binding site of the RBD and hACE2 is competitive. Several studies have been conducted using competitive inhibitors against the main protease from SARS-CoV-2 [77]. This competitive inhibition of bromelain-derived peptide can prevent the adhesion of hACE2 to the RBD since the positions of key amino acids in the RBD have been filled by the bromelain-derived peptide. This was confirmed by He et al. [78], who stated that a competitive inhibitor that binds to the active site of an enzyme can inhibit the activity of the enzyme by competing with the substrate, which in this case is hACE2. Due to the similarity of their binding interfaces to the hACE2 receptor's S1 binding site (PDB ID: 6M0J), the antimicrobial peptide dermaseptin and its analog have been shown to act as competitive inhibitors for the hACE2 receptor [79].

Molecular Dynamics Simulations Study
The beneficial uses of MD simulation include developing a greater understanding of protein-ligand interactions, determining the spatial orientation of active sites, determining the motion of active site residues, and analyzing protein conformational dynamics. MD simulations with all atoms were run for 100 ns on the peptide-protein complexes of bromelain-derived peptide and the RBD. To understand the deviation of Cα atoms from the backbone as well as the fluctuations in amino acid contact during simulation, RMSD and RMSF analyses were conducted.
Currently, MD simulations are widely used in drug and vaccine design studies [80,81]. The notion that in silico techniques will aid in the treatment of diseases in desperate need of a cure, such as COVID-19, is becoming increasingly important. In particular, techniques such as molecular docking and MD simulation can save time and costs in demanding drug development efforts [82]. MD simulations are frequently used to investigate the stability of protein-ligand, protein-peptide, protein-protein, and protein-DNA/RNA complexes [83]. By MD simulation, the stability of the interactions between SARS-CoV-2 S1b unit (RBD) variants and the bromelain-derived peptide was investigated by employing the HADDOCK server. Amber99sb-ildn was chosen since it is suitable for force field protein-peptide and protein-protein simulations. A total of 10 protein-peptide complexes obtained by the molecular docking method were subjected to a pre-simulation for 20 ns. In this presimulation, WT and RBD variants with high stability were investigated over a longer time via a 100 ns simulation.
The stability of the complex of bromelain peptide and RBDs was demonstrated by RMSD, Rg, and RMSF trajectory analysis. RMSD measurement is the main parameter expressing the deviation and shift in protein structure. The RMSD value of the complexes was calculated according to the backbone atoms. Rg calculations are another important parameter that provides information about the compactness of a protein. The lower and more constant the RMSD and Rg values are, the more stable the complex structure is. As shown in Figure 6, the protein-peptide complexes remained stable after a certain period of time. The RMSD value of the BR/RBD-peptide complex was approximately 0.3 nm, and the Rg value remained constant between 1.76 and 1.86 nm. The RMSD value of the US/RBD-peptide complex slowly increased up to 60 ns and approached 0.3 nm, remaining stable thereafter. The Rg value varied between 1.82 nm and 1.88 nm throughout the simulation. In the WT/RBD-peptide complex, it rose slightly above 0.3 ns during the first 10 ns and then remained stable at approximately 0.2 nm after 20 ns. Similar to the RMSD chart, the Rg chart shows fluctuation and stability. The highest Rg value was measured at 1.90 nm. The RMSD value of the UK/RBD-peptide complex increased to 0.3 nm in the first 20 ns and then remained constant at approximately 0.2 nm until 80 ns, after which the protein-peptide interaction was lost. Rg remained constant between 1.85 and 1.90 nm until 80 ns and increased greatly thereafter. Finally, the NG/RBD-peptide assembly remained stable below 0.2 nm throughout the simulation, with the Rg value also measuring between 1.85 and 1.87 nm.
RMSF is another important MD trajectory analysis parameter showing fluctuations and conformational changes in protein structure. In the residue-based RMSF analysis, the numerical values of amino acids with high mobility are high, while the RMSF values of residues with lower mobility are low. RMSF calculations were made based on protein C-α atoms. To evaluate the status of the complex formed by bromelain peptide with RBD WT and RBD variants, the peptide-free WT/RBD (also known as the apo form) was simulated under the same conditions. As seen in Figure 7, the apo form has more mobility when compared to the amino acid holo forms (numbered 346, 382-386, 389, 454-456, and 517-519). Since the variants and WT/RBD-peptide complexes became more stable than the apo form, the peptide-protein complexes can be said to remain stable. The RBD flexibility of the BR/RBD (0.30 nm) and UK/RBD (0.21 nm) variants, which contain the N501Y mutation, was higher than that of WT/RBD (0.09 nm). This may explain why the infectiousness of the N501Y mutant is higher. The E484K mutation also increases the mobility of the RBD. In particular, the flexibility of BR/RBD (0.24 nm) and NG/RBD (0.21 nm) variants containing E484K increased at amino acid 484 compared to WT/RBD (0.13 nm).

MM-PBSA Binding-Free Energy
Although BFE calculations were performed in the preceding section, we consider calculations based on MM-PBSA in this section. The MM-PBSA was calculated to measure the in-depth atomic-level interaction energy of bromelain-derived peptide-RBD complexes ( Table 4). This is due to the lower accuracy of the BFE value derived from docking results [84], which is only used as an initial prediction. It employs a scoring function to rank the various possible poses of a ligand in a binding pocket, which is focused on determining the binding affinity. The BR and NG variants with the most stable RMSD plots and the bromelain-derived peptide produced the lowest average binding energy to the peptide (−287.356, and −255.801 kJ/mol, respectively). Similarly, according to the RMSD plot, the UK variant, with the lowest stability and disrupted protein-peptide complex after 80 ns, gave an average free energy value of −89.129 kJ/mol. The WT and US variants gave values of −173.243 and −150.460 kJ/mol, respectively. According to these MM-PBSA results, the protein-peptide connection continues for a certain period of time. This may indicate that the bromelain peptide remains stable at the possible binding site to the RBD.

Conclusions
In this study, we evaluated bromelain-derived peptide as a potential drug to target inhibition of the interaction between the RBD and the hACE2 receptor. The results of the MD simulation validated our decision to propose the bromelain-derived peptide as an inhibitor against the RBD WT and its variants. Throughout the simulation, the peptide-RBD complexes remained stable. This result indicates that bromelain-derived peptide could potentially be used as a drug to prevent a prolonged COVID-19 pandemic by inhibiting viral fusion and entry into cells. However, further in vitro and in vivo testing will be required to validate the efficacy and safety of the bromelain-derived peptide as a SARS-CoV-2 inhibitor.
Supplementary Materials: The following are available online, Figure S1: Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelain-derived peptide and RBD SA by LigPlot+. RBD SA contains the mutations K417N, E484 and N501Y. Figure S2: Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelain-derived peptide and RBD UK by LigPlot+. RBD UK contains the mutation N501Y. Figure S3: Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelain-derived peptide and RBD CA (US) by LigPlot+. RBD CA contains the mutation L452R. Figure S4. Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelainderived peptide and RBD SG (NY1) by LigPlot+. RBD SG contains the mutation S477G. Figure S5. Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelain-derived peptide and RBD SN (NY2) by LigPlot+. RBD SN contains the mutation S477N. Figure S6. Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelain-derived peptide and RBD SC by LigPlot+. RBD SC contains the mutation N439K. Figure S7. Sketch diagram depicting the 3D (A) and 2D (B) interactions between bromelain-derived peptide and RBD NG by LigPlot+. RBD NG contains the mutation E484K. Figure  Data Availability Statement: The authors will make the raw data supporting the conclusions of this manuscript available to any qualified researcher.